function StrainZero(GSS)
%StrainZero - Strain computation in selected points
% 
%       StrainZero(GSS)
%
% Function for the calculation of strain in selected points.
% The input argument is a GSS variable provided in a GridStrain session.
% No more than data about positions, velocities, velocity erros (and 
% error ellipses) and (optionally) excluded points are required in GSS 
% (all other data in GSS, if available, are ignored by StrainZero).
% If GSS is undefined or empty, the name of a file with such a variable
% can be interactively managed.
% 
% The computations can be carried out either on the experimental points
% (i.e. the stations whose displacements or velocities are known), or  
% on interactively defined points.
% If the computation is carried out on the experimental points, the 
% chi-square and the reduced chi-square are also computed and shown.
% If one or more experimental points were excluded from the computations
% in the GridStrain session which provided the used GSS variable, the 
% chi-square and the reduced chi-square can be computed eiter for either 
% all stations or kept stations only. 
%
% For each point in which the computation can be carried out, the following
% output data are provided on the screen:
%
%   POINT COORDINATES;
%   MAXIMUM AND MINIMUM STRAINS (OR STRAIN RATES);
%   ERRORS ON STRAINS (OR STRAIN RATES);
%   ANGLE OF THE MAXIMUM STRAIN DIRECTION WITH RESPECT TO THE NORTH
%       DIRECTION;
%   CHANGE IN AREA (OR RATE OF CHANGE IN AREA);
%   NORMALIZED SHEAR (OR RATE OF NORMALIZED SHEAR);
%   GEOMETRIC SIGNIFICANCE OF THE RESULTS;
%   MODELED TRANSLATIONS (OR VELOCITIES);
%   ERRORS ON MODELED TRANSLATIONS (OR VELOCITIES); 
%   RIGID BODY ROTATION AND CORRESPONDING ERROR  
%
% The results can be saved on an ASCII file.

% A. Pesci, G. Teza, 2006, 2008, 2011, 2022.

Options = GridStrainOptions;

if nargin < 1
    
    [filen, pathn] = uigetfile({'*.mat','MATLAB (*.mat) files'},...
        'INPUT MATLAB FILE');
    if isequal(filen,0) || isequal(pathn,0)
        fprintf('\nNo a data file was selected \n');
        return
    else
        filena = fullfile(pathn,filen);
    end
    GSSLoad = load(filena);
    try
        GSS = GSSLoad.GSS;
        if isempty(GSS.XGridStep)
            Iok = false;
        else
            Iok = true;
        end
    catch
        Iok = false;
    end
    if ~Iok
        fprintf('\nInvalid input data \n');
        return
    end
   
end

ns = GSS.StationID;
stn = GSS.StationName;
e = GSS.East;
n = GSS.North;
ve = GSS.EastVelocity;
vn = GSS.NorthVelocity;
eve = GSS.EastVelocityError;
evn = GSS.NorthVelocityError;
a = GSS.ErrorEllipsea;
b = GSS.ErrorEllipseb;
theta = GSS.ErrorEllipsetheta;
IMS = GSS.RefImageStruct;
try 
    npe = GSS.IncludedPoints;
catch
    npe = true(size(e));
end
if isempty(ns)
    nso = (1:numel(e))';
else
    nso = ns;
end

mencalc = menu('CALCULATION OPTIONS:',...
    'USER-DEFINED SCALE FACTOR',...
    'INFINITE SCALE FACTOR (Uniform strain model)');
if mencalc == 1
    d0 = mean([max(e)-min(e),max(n)-min(n)])/5;
    prompty = {'Scale factor'};
    defansy = {num2str(d0)};
    fields = {'d0'};
    acquis = inputdlg(prompty,'SCALE FACTOR MANAGEMENT',1,defansy);
    if ~isempty(acquis)
        acq = cell2struct(acquis,fields);
        d0m = str2double(acq.d0);
    else
        d0m = [];
    end
    if ~isempty(d0m)
        d0 = d0m; 
    else
        d0 = Inf; 
    end
else
    d0 = Inf;
end

Ii = npe;       % npe = 1 <-> point included
le = length(e);
Ix = ~npe;
ei  = e; 
vei = ve; 
ex  = e; 
vex = ve;
ni  = n; 
vni = vn; 
nx  = n; 
vnx = vn;
ei(Ix)  = NaN; 
vei(Ix) = NaN; 
ex(Ii)  = NaN; 
vex(Ii) = NaN;
ni(Ix)  = NaN; 
vni(Ix) = NaN; 
nx(Ii)  = NaN; 
vnx(Ii) = NaN;
eve(Ix) = eve(Ix)*100000;
evn(Ix) = evn(Ix)*100000;  
if ~isempty(ns)
    nsi = ns;
    nsi(Ix) = NaN;
    nsx = ns;
    nsx(Ii) = NaN;
    stni = [];
    stnx = [];
else
    nsi = [];
    nsx = [];
    stni = stn;
    stni(Ix) = {[]};
    stnx = stn;
    stnx(Ii) = {[]};
end

figs = figure; 
set(figs,'Units','normalized','OuterPosition',[0 0 1 1]); 
plot(ei,ni,'ok',ex,nx,'og');
if Options.VisualizePoints
    plot(ei,ni,'ok',ex,nx,'og'); 
    if Options.VisualizePointIDs
        addIDName(ei,ni,vei,vni,nsi,stni,ex,nx,vex,vnx,nsx,stnx);
    end
end
legend('Point included','Point not included','Location','NorthEastOutside');
hold on;
quiver(ei,ni,vei,vni,'k','HandleVisibility','off'); 
quiver(ex,nx,vex,vnx,'m','HandleVisibility','off');
title('VELOCITY (OR DISPLACEMENT) FIELD AND STRAIN','FontSize',14);
xlabel('EAST','FontSize',14); 
ylabel('NORTH','FontSize',14); 
axis equal; 
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',12);
mine = min(e); 
maxe = max(e);
minn = min(n); 
maxn = max(n);

if sum(Ii) < le  % case of excluded points (Ii is the set of indices of included points, see above)
    menchof = menu('STRAIN RATE COMPUTATION GENERAL OPTION:',...
        'COMPUTATION ON ALL THE EXPERIMENTAL POINTS',...
        'COMPUTATION ON THE KEPT EXPERIMENTAL POINTS ONLY',...
        'COMPUTATION ON INTERACTIVELY SELECTED POINTS');
else                % case of all kept points  
    menchof = menu('STRAIN RATE COMPUTATION GENERAL OPTION:',...
        'COMPUTATION ON THE EXPERIMENTAL POINTS',...
        'COMPUTATION ON INTERACTIVELY SELECTED POINTS');
    if menchof == 2
        menchof = 3; 
    end
end

if ismember(menchof,[1 2]) 
    mout = zeros(le,21);
    mout(:,1:10) = [nso e n ve vn eve evn a b theta];
    ec = e; 
    nc = n; 
    lec = le;
    if menchof == 2 
        mout = mout(Ii,:);
        if ~isempty(stn)
            stn = stn(Ii);
        end
        ec = e(Ii); 
        nc = n(Ii);
        lec = numel(ec);
    end
else
    mout = [];
end

nopo = 0;
kc   = 1;  % index of analyzed point
kcvv = 1;  % index for the scale arrow visualization
while nopo == 0
    
    if ismember(menchof,[1 2])
        menpo = 4;   % case initial GPS stations
    else
        menpo = menu('CHOICE OF ESTIMATION POINT:',...
            '1 - BY COORDINATES',...
            '2 - BY CLICK ON PLOT',...
            '3 - CENTER OF MASS');
    end

    if menpo == 1   % by coordinates
        
        xpd = mean([mine maxe]);
        ypd = mean([minn maxn]);
        prompty = {'x coordinate:','y coordinate'};
        defansy = {num2str(xpd),num2str(ypd)};
        fields  = {'xp','yp'};
        acquis  = inputdlg(prompty,'CHOICE OF POINT',1,defansy);
        if ~isempty(acquis)
            acq = cell2struct(acquis,fields);
            xp  = str2double(acq.xp);
            yp  = str2double(acq.yp);
            if (xp < mine)||(xp > maxe)||(yp < minn)||(yp > maxn)
                menwa = menu(...
                    'WARNING: SELECTED POINT IS OUTSIDE THE DATA ZONE',...
                    'CONTINUE CONSIDERING THIS POINT',...
                    'CHOICE OF ANOTHER POINT',...
                    'EXIT FROM STRAINZERO');
            end
        else
            xp = xpd;
            yp = ypd;
            menwa = 1;  % continue with this point 
        end
        
    elseif menpo == 2   % by click on plot 
        
        figure(figs);
        [xp,yp] = ginput(1);
        menwa = 1;  % For the continuation of computation
    
    elseif menpo == 3   % center of mass
        
        we = (1./eve).^2; 
        wn = (1./evn).^2;
        xp = dot(we,e)/sum(we);
        yp = dot(wn,n)/sum(wn);
        menwa = 1;
    
    else                % experimental points
        
        if kc <= lec
            xp = ec(kc);
            yp = nc(kc);
            menwa = 1;
        else
            menwa = 3;
        end
        
    end
    
    if menwa == 1
        
        %------------------------------------------------------------------
        [Emax,Emin,Phi,eEmax,eEmin,n0,flc,fls,u0,v0,eu0,ev0,...
            omeg,eomeg] = gs_strain(xp,yp,e,n,ve,vn,eve,evn,...
            npe,d0,Options); %#ok<ASGLU>
        %------------------------------------------------------------------
        
        xfv = [NaN; xp]; % NaN is necessary for the correct visualization  
        
        if flc == 1
            
            yfv = [NaN; yp];
            if kcvv == 1
                mama  = max(abs(Emax));
                mami  = max(abs(Emin));
                diffe = max(e)-min(e);
                diffn = max(n)-min(n);
                scalfact = 0.25*min([diffe diffn])/max(mama,mami);
            end
            
            if Phi >= 0
                phi = Phi;
            else
                phi = pi+Phi;
            end
            
            EIG1B = nan(2); % In order to draw the first eigenvetor
            EIG1R = nan(2); % B (blue): dilatation; R (red): compression   
            EIG2B = nan(2);
            EIG2R = nan(2);
            aema  = abs(Emax);
            aemi  = abs(Emin);
            cphi  = cos(phi);
            sphi  = sin(phi);
            
            if phi <= pi/2   % Angle in the range 0-90�
                
                if Emax >= 0
                    EIG1B(2,1) = aema*sphi;
                    EIG1B(2,2) = aema*cphi;
                else
                    EIG1R(2,1) = aema*sphi;
                    EIG1R(2,2) = aema*cphi;
                end
                if Emin >= 0
                    EIG2B(2,1) = aemi*cphi;
                    EIG2B(2,2) = -aemi*sphi;
                else
                    EIG2R(2,1) = aemi*cphi;
                    EIG2R(2,2) = -aemi*sphi;
                end
                
            else
                
                if Emax >= 0
                    EIG1B(2,1) = aema*sphi;  % sin(pi-phi) = sin(phi)
                    EIG1B(2,2) = aema*cphi;  % cos(pi-phi) = -cos(phi) 
                else
                    EIG1R(2,1) = aema*sphi;
                    EIG1R(2,2) = aema*cphi;
                end
                if Emin >= 0
                    EIG2B(2,1) = -aemi*cphi; 
                    EIG2B(2,2) = aemi*sphi;
                else
                    EIG2R(2,1) = -aemi*cphi;
                    EIG2R(2,2) = aemi*sphi;
                end
            
            end
           
            figure(figs);
            hold on;
            quiver(xfv,yfv,...
                scalfact*EIG1B(:,1),scalfact*EIG1B(:,2),...
                0,'b','HandleVisibility','off');
            quiver(xfv,yfv,...
                -scalfact*EIG1B(:,1),-scalfact*EIG1B(:,2),...
                0,'b','HandleVisibility','off');
            quiver(xfv-scalfact*EIG1R(:,1),yfv-scalfact*EIG1R(:,2),...
                scalfact*EIG1R(:,1),scalfact*EIG1R(:,2),...
                0,'r','HandleVisibility','off');
            quiver(xfv+scalfact*EIG1R(:,1),yfv+scalfact*EIG1R(:,2),...
                -scalfact*EIG1R(:,1),-scalfact*EIG1R(:,2),...
                0,'r','HandleVisibility','off');
            quiver(xfv,yfv,...
                scalfact*EIG2B(:,1),scalfact*EIG2B(:,2),...
                0,'b','HandleVisibility','off');
            quiver(xfv,yfv,...
                -scalfact*EIG2B(:,1),-scalfact*EIG2B(:,2),...
                0,'b','HandleVisibility','off');
            quiver(xfv-scalfact*EIG2R(:,1),yfv-scalfact*EIG2R(:,2),...
                scalfact*EIG2R(:,1),scalfact*EIG2R(:,2),...
                0,'r','HandleVisibility','off');
            quiver(xfv+scalfact*EIG2R(:,1),yfv+scalfact*EIG2R(:,2),...
                -scalfact*EIG2R(:,1),-scalfact*EIG2R(:,2),...
                0,'r','HandleVisibility','off');
            
            if kcvv == 1    % units representation: 
                amax  = abs(Emax);
                mamax = sprintf('%0.0g',amax);
                quiver(maxe+0.2*diffe,minn,scalfact*amax,0,0,'b',...
                    'HandleVisibility','off');
                text(maxe+0.1*diffe,minn+0.025*diffn,['STRAIN SCALE: ' mamax]);
            end
            kcvv = kcvv+1;
            axis equal; 
            hold off;
            
            %--------------------------------------------------------------
            rep_point = gs_reppoint(xp,yp,Emax,Emin,eEmax,eEmin,Phi,d0,...
                fls,u0,v0,eu0,ev0,omeg,eomeg);
            %--------------------------------------------------------------
            
        else
            rep_point = sprintf('POINT: \n East: %f, \n North: %f;\n %s.',...
                xp,yp,'=== IN THIS POINT, CALCULATIONS CANNOT BE PERFORMED ===');
        end
        fprintf('\n-----------------------------------------------------------');
        disp(rep_point);
        fprintf('-----------------------------------------------------------\n');
        
        if ismember(menchof,[1 2])
            
            mout(kc,11) = Emax;  mout(kc,12) = Emin;
            mout(kc,13) = eEmax; mout(kc,14) = eEmin;
            mout(kc,15) = Phi;
            mout(kc,16) = u0;    mout(kc,17) = v0;
            mout(kc,18) = eu0;   mout(kc,19) = ev0;
            mout(kc,20) = omeg;  mout(kc,21) = eomeg;
            
            if kc >= lec
                nopo = 1; 
            end
        
        else
            
            mout(kc,1)  = kc;
            mout(kc,2)  = xp;    mout(kc,3)  = yp;
            mout(kc,4)  = Emax;  mout(kc,5)  = Emin;
            mout(kc,6)  = eEmax; mout(kc,7)  = eEmin;
            mout(kc,8)  = Phi;
            mout(kc,9)  = u0;    mout(kc,10) = v0;
            mout(kc,11) = eu0;   mout(kc,12) = ev0;
            mout(kc,13) = omeg;  mout(kc,14) = eomeg;
            
            menuc = menu(rep_point,'SELECT ANOTHER POINT','SKIP');
            if menuc == 2
                nopo = 1; 
            end
        
        end
        
    elseif menwa == 3       % exit from StrainZero (see above definition of menwa)
        nopo = 1;
    end
    
    kc = kc+1;

end

if ismember(menchof,[1 2])   % cases where chi-square should be computed 
    
    if menchof == 1
        Ii = 1:le; 
        lec = le; 
    end
    xs = zeros(2*lec,1); 
    exs = zeros(2*lec,1); 
    xm = zeros(2*lec,1);
    vec = ve(Ii); 
    vnc = vn(Ii);
    xs(1:lec)  = vec;  
    xs(lec+1:end) = vnc;
    exs(1:lec) = eve(Ii); 
    exs(lec+1:end) = evn(Ii);
    xm(1:lec)  = mout(:,16); 
    xm(lec+1:end) = mout(:,17);
    [chi2,chi2r,dof] = chi2calc(xs,xm,exs,6);
    fprintf(...
        '\nCHI-SQUARE COMPUTATION (SCALE FACTOR: %f):\n COMPLETE: %5.2f;\n REDUCED: %5.2f; \n DOFs: %d\n\n',...
        d0,chi2,chi2r,dof);
    arry_vis(ec,nc,vec,vnc,mout(:,16),mout(:,17),ns,stn,...
        'INITIAL (BLUE) AND MODELED (RED) VECTORS',IMS,Options);
 
elseif ~isempty(mout)
    
    arry_vis(mout(:,2),mout(:,3),mout(:,9),mout(:,10),[],[],ns,stn,...
        'MODELED DISPLACEMENT OR VELOCITY VECTORS',IMS,Options);
    
end
    
if ~isempty(mout)
    mensav = menu('DATA SAVING','YES','NO');
    if mensav == 1
        [fileno, pathno] = uiputfile(...
            '*.txt','SAVE STRAIN_ZERO RESULTS (ASCII FILE)','untitled.txt');
        filenout = fullfile(pathno,fileno);
        if isempty(stn)
            save(filenout,'mout','-ascii','-double');
        else
            moutr = mout(:,2:end);
            [nrout,ncout] = size(mout);
            % string for fprintf (note that %f to be added are ncout -1):
            sout = '%s\t'; 
            for kk = 1:(ncout-1)
                if kk < (ncout-1)
                    sout = [sout '%f\t']; %#ok<*AGROW>
                else
                    sout = [sout '%f\n'];
                end
            end
            fid = fopen(filenout,'w+');
            for k = 1:nrout
                fprintf(fid,sout,stn{k},moutr(k,:));
            end
            fclose(fid);
        end
    end
end


%% ANCILLARY FUNCTIONS 

function addIDName(ei,ni,vei,vni,nsi,stni,ex,nx,vex,vnx,nsx,stnx) 
% Add station ID or station name
mediff = mean([nanstd(diff(ei)) nanstd(diff(ni))]); %#ok<NANSTD>
for kk = 1:numel(ei)
    if ~isempty(nsi) && ~isnan(nsi(kk))
        strikk = num2str(nsi(kk));
        okikk = true;
    elseif ~isempty(stni) && ~isempty(stni{kk})
        strikk = stni{kk};
        okikk = true;
    else
        okikk = false;
    end
    if okikk 
        se = sign(vei(kk)); 
        sn = sign(vni(kk));
        ht = text(ei(kk)-se*mediff/12,ni(kk)-sn*mediff/12,strikk);
        set(ht,'Color','k');
    end
end
for kk = 1:numel(ex)
    if ~isempty(nsx) && ~isnan(nsx(kk))
        strxkk = num2str(nsx(kk));
        okxkk = true;
    elseif ~isempty(stnx) && ~isempty(stnx{kk})
        strxkk = stnx{kk};
        okxkk = true;
    else
        okxkk = false;
    end
    if okxkk 
        se = sign(vex(kk)); sn = sign(vnx(kk));
        ht = text(ex(kk)-se*mediff/12,nx(kk)-sn*mediff/12,strxkk);
        set(ht,'Color','m');
    end
end

function arry_vis(e,n,ve,vn,ve1,vn1,ns,stn,str_title,IMS,Options)
% Quiver and arrow scale visualization on the points (e,n) with the
% velocities ve, vn, shown as blue arrows, and the optional velocities
% ve1, vn1 (if undefined, they must be empty vectors), shown as red arrows.
e = e(:); 
n = n(:); 
ve = ve(:); 
vn = vn(:); 
ve1 = ve1(:); 
vn1 = vn1(:);
% A new point is generated in order to represent the scale arrow
esca = e;  
nsca = n; 
lv = length(e);  
vesca = ve; 
vnsca = vn;                 
std_esca = std(esca-mean(esca));        
std_nsca = std(nsca-mean(nsca));
esca(lv+1)  = min(esca)+std_esca/10;
nsca(lv+1)  = min(nsca)+std_nsca/10;
if isempty(ve1)
    vesca(lv+1) = (max(abs(ve))+max(abs(vn)))/4;
    scalef = 0.2*max(max(e)-min(e),max(n)-min(n))/max(max(abs(ve),abs(vn)));
else
    ve2 = [ve; ve1]; vn2 = [vn; vn1];
    vesca(lv+1) = (max(abs(ve2))+max(abs(vn2)))/4;
    scalef = 0.2*max(max(e)-min(e),max(n)-min(n))/max(max(abs(ve2),abs(vn2)));
end
valnum = vesca(end);
valnum = valnum*1000;               % to always have data in mm
val = fixed_dig(num2str(valnum),4); % scale vector intensity
vnsca(lv+1) = 0.00;                 % horizontal vector
figs = figure; 
set(figs,'Units','normalized','OuterPosition',[0 0 1 1]); 
plot(e,n,'xr');
quiver(esca,nsca,vesca*scalef,vnsca*scalef,0,'b');
if Options.VisualizePointIDs
    addIDName(e,n,ve,vn,ns,stn,[],[],[],[],[],[]);
end
if ~isempty(ve1)
    hold on;
    quiver(e,n,ve1*scalef,vn1*scalef,0,'r');
    hold off;
end
text(esca(end),nsca(end)+std_nsca/10,[val ' mm']);
title(str_title,'FontSize',14);
xlabel('EAST','FontSize',14); 
ylabel('NORTH','FontSize',14); 
axis equal;
ALA = axis;
gs_ReadBGImage(IMS,ALA);
hold off; 
set(gca,'FontSize',12);